Aim: Demonstrate photometry on a series of bias and flat field corrected images of a Near Earth Asteroid.
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
## make matplotlib appear in the notebook rather than in a new window
%matplotlib inline
In [ ]:
datadir = ''
objname = '2016HO3'
In [ ]:
def plotfits(imno):
img = fits.open(datadir+objname+'_{0:02d}.fits'.format(numb))[0].data
f = plt.figure(figsize=(10,12))
#im = plt.imshow(img, cmap='hot')
im = plt.imshow(img[480:580, 460:600], cmap='hot')
plt.clim(1800, 2800)
plt.colorbar(im, fraction=0.034, pad=0.04)
plt.savefig("figure{0}.png".format(imno))
plt.show()
In [ ]:
numb = 1
plotfits(numb)
In [ ]:
numb = 2
plotfits(numb)
Select part of the image for ease of display.
In [ ]:
partimg = fits.open(datadir+objname+'_01.fits'.format(numb))[0].data[480:580, 460:600]
Define starting values. Fill in values here:
In [ ]:
targcen = np.array([##,##]) ## target center
compcen = np.array([##,##]) ## comparison center
Aperture photometry set up. Play around with adjusting the aperture radii sizes and see the resulting image under 'Tests'
In [ ]:
searchr = 6 ## search box size
ap_r = 2 ## aperture radius
sky_inner = 3
sky_outer = 5
Calculate Center of Mass (CoM) defined as: $\bar{x} = \frac{\sum A_i x_i}{\sum A_i }$, $\bar{y} = \frac{\sum A_i y_i}{\sum A_i }$.
In [ ]:
def cent_weight(n):
"""
Assigns centroid weights
"""
wghts=np.zeros((n),np.float)
for i in range(n):
wghts[i]=float(i-n/2)+0.5
return wghts
def calc_CoM(psf, weights):
"""
Finds Center of Mass of image
"""
cent=np.zeros((2),np.float)
### Write Equations for finding Center of Mass here ###
return cent
Use centroiding algorithm to find the actual centers of the targe and comparison.
In [ ]:
## Cut a box between search limits, centered around targcen
targbox = partimg[targcen[0]-searchr : targcen[0]+searchr, targcen[1]-searchr : targcen[1]+searchr]
weights = cent_weight(len(targbox))
tcenoffset = calc_CoM(targbox, weights)
print(tcenoffset)
tcenter = targcen + tcenoffset
Inspect PSF to see whether shift makes sense
In [ ]:
plt.plot(sum(targbox))
plt.show()
In [ ]:
compbox = partimg[compcen[0]-searchr : compcen[0]+searchr, compcen[1]-searchr : compcen[1]+searchr]
compw = cent_weight(len(compbox))
ccenoffset = calc_CoM(compbox,compw)
ccenter = compcen + ccenoffset
In [ ]:
print(tcenter)
In [ ]:
compw
In [ ]:
def circle(npix, r1):
"""
Builds a circle
"""
pup=np.zeros((npix,npix),np.int)
for i in range(npix):
for j in range(npix):
r=np.sqrt((float(i-npix/2)+0.5)**2+(float(j-npix/2)+0.5)**2)
if r<=r1:
pup[i,j]=1
return pup
In [ ]:
def annulus(npix, r_inner,r_outer=-1.):
"""
Builds an annulus
"""
pup=np.zeros((npix,npix),np.int)
for i in range(npix):
for j in range(npix):
#### Fill in annulus form here ####
if ((r<=r_outer)&(r>=r_inner)):
pup[i,j]=1
return pup
Create mask
In [ ]:
circmask = circle(searchr*2, ap_r)
annmask = annulus(searchr*2, sky_inner, sky_outer)
Define new regions where the target and comparison are centered.
In [ ]:
newtarg = partimg[int(round(tcenter[0]))-searchr : int(round(tcenter[0]))+searchr, int(round(tcenter[1]))-searchr : int(round(tcenter[1]))+searchr]
newcomp = partimg[int(round(ccenter[0]))-searchr : int(round(ccenter[0]))+searchr, int(round(ccenter[1]))-searchr : int(round(ccenter[1]))+searchr]
Place mask on region
In [ ]:
targaper = newtarg * circmask
compaper = newcomp * circmask
Place mask on sky annulus slice.
In [ ]:
targann = newtarg * annmask
compann = newcomp * annmask
a. Display image with target and comparison centers before and after centroiding
In [ ]:
im = plt.imshow(partimg, cmap='hot')
plt.clim(1800, 2800)
plt.scatter(targcen[1], targcen[0], c='g', marker='x')
plt.scatter(compcen[1], compcen[0], c='g', marker='x')
plt.scatter(tcenter[1], tcenter[0], c='b', marker='x')
plt.scatter(ccenter[1], ccenter[0], c='b', marker='x')
plt.show()
b. Disply image with aperture mask and sky annulus
In [ ]:
im = plt.imshow(targaper, cmap='hot')
plt.clim(1800, 2800)
plt.show()
In [ ]:
im = plt.imshow(targann, cmap='hot')
plt.clim(1800, 2800)
plt.show()
Calculate Signal-to-Noise Ratio. CCD noise = sqrt(signal + background + dark current + read noise). Ignore dark current and read noise here.
In [ ]:
def calcsnr(target, bg):
signal = target - bg
noise = np.sqrt(signal + bg)
snr = signal / noise
return snr, noise
Sum all flux inside target and comparison apertures and divide by number of pixels to get average count per pixel.
In [ ]:
targc = np.sum(targaper) / np.sum(circmask)
targbg= np.sum(targann) / np.sum(annmask)
compc = np.sum(compaper) / np.sum(circmask)
compbg= np.sum(compann) / np.sum(annmask)
In [ ]:
snr, noise = calcsnr(targc, targbg)
print(snr)
In [ ]:
snr, noise = calcsnr(compc, compbg)
print(snr)
In [ ]:
## Write code here that tries a range of photometry apertures and finds the best SNR ##
In [ ]:
print(bestaper)
print(snr)
Given the comparison is of known magnitude of 19.4
In [ ]:
targc = circle(searchr*2, ap_r)*newtarg
targskyc = annulus(searchr*2, sky_inner, sky_outer)*newtarg
compc = circle(searchr*2, ap_r)*newcomp
compskyc = annulus(searchr*2, sky_inner, sky_outer)*newcomp
ratio = np.sum(compc)/np.sum(targc)
### complete here ###
### complete here ###
### complete here ###
refmag = 19.4
### complete here ###
print("Measured Magnitude = {:0.3f} ± {:0.3f}".format(mag, sigmamag))